pw = @(lon) (lon-min([dc.lon]))/(max([dc.lon])-min([dc.lon]))*(1-2*m)+m;
ph = @(lat) (lat-min([dc.lat]))/(max([dc.lat])-min([dc.lat]))*(1-2*m)+m;
figure('position',[0 0 1000 1000]);
plot(coastlon,coastlat,'-k');
xlim([min([dc.lon]) max([dc.lon])])
ylim([min([dc.lat]) max([dc.lat])])
[esp,hsp,~,~,Mdd(i_d,1)] = histvdens(dc(i_d).dd(id_sp,:), dens(id_sp,:,i_d).*v_g(id_sp,:,i_d), 100);
[eau,hau,~,~,Mdd(i_d,2)] = histvdens(dc(i_d).dd(id_au,:), dens(id_au,:,i_d).*v_g(id_au,:,i_d), 100);
%[e2016,h2016,~,~,Mdd(i_d,3)] = histvdens(dc(i_d).dd(id_2016,:), dens(id_2016,:,i_d).*v_g(id_2016,:,i_d), 100);
subAX = axes('position',[pw(dc(i_d).lon)-w/2 ph(dc(i_d).lat)-h/2 w h]);
%polarplot(deg2rad(e2016),h2016);hold on
polarplot(deg2rad(esp),hsp,'Color',id_sp_col,'linewidth',2);hold on
polarplot(deg2rad(eau),hau,'Color',id_au_col,'linewidth',2);
%tmp=round(reshape(dens(id_sp,:,i_d).*v_g(id_sp,:,i_d),1,[])); tmp(isnan(tmp)) = 0;
%polarhistogram(repelem(deg2rad(reshape(dc(i_d).dd(id_sp,:),1,[])), tmp),38); hold on;
%polarplot(angle([Mdd(i_d,1) Mdd(i_d,1)]),[0 max(hsp)],'--','Color',id_sp_col,'LineWidth',2)
%polarplot(angle([Mdd(i_d,2) Mdd(i_d,2)]),[0 max(hau)],'--','Color',id_au_col,'LineWidth',2)
% polarplot(angle([Mdd(i_d,3) Mdd(i_d,3)]),[0 max([h2016;h2016])],'LineWidth',2)
Ax = gca; Ax.ThetaGrid = 'off'; Ax.RGrid='off'; Ax.RTickLabel = []; Ax.ThetaTickLabel = [];
title(dc(i_d).name,'Color',id_sno_col(id_no(i_d)+1,:))
figure('position',[0 0 1200 1000]);
subplot(2,3,[1 2]); hold on;
[b,h,Mg,Sg] = histvdens(v_g(id_sp,:,:),dens(id_sp,:,:));
bar(b,h,'FaceAlpha',0.6, 'DisplayName',['Groundspeed. M=' num2str(Mg,3) ' ; S=' num2str(Sg,3)]);
[b,h,Ma,Sa] = histvdens(v_a(id_sp,:,:),dens(id_sp,:,:));
bar(b,h,'FaceAlpha',0.6, 'DisplayName',['Airspeed. M=' num2str(Ma,3) ' ; S=' num2str(Sa,3)]);
legend; xlabel(['Spring. Bird speed [m/s]. Total nb: ' num2str(nansum(nansum(nansum(dens(id_sp,:,:))))/1000000,3) ' M']);
ylabel('Sum of bird density [bird/km^3]'); xlim([0 30]); box on
subplot(2,3,[4 5]); hold on;
[b,h,Mg,Sg] = histvdens(v_g(id_au,:,:),dens(id_au,:,:));
bar(b,h,'FaceAlpha',0.6, 'DisplayName',['Groundspeed. M=' num2str(Mg,3) ' ; S=' num2str(Sg,3)]);
[b,h,Ma,Sa] = histvdens(v_a(id_au,:,:),dens(id_au,:,:));
bar(b,h,'FaceAlpha',0.6, 'DisplayName',['Airspeed. M=' num2str(Ma,3) ' ; S=' num2str(Sa,3)]);
legend; xlabel(['Autumn. Bird speed [m/s]. Total nb: ' num2str(nansum(nansum(nansum(dens(id_au,:,:))))/1000000,3) ' M']);
ylabel('Sum of bird density [bird/km^3]'); xlim([0 30]); box on
subplot(4,3,(i_s-1)*3+3); hold on
[b,h,Mg,Sg] = histvdens(v_g(id_s(i_s,:),:,:),dens(id_s(i_s,:),:,:));
bar(b,h,'FaceAlpha',0.6, 'DisplayName',['M=' num2str(Mg,3) ' ; S=' num2str(Sg,3)]);
[b,h,Ma,Sa] = histvdens(v_a(id_s(i_s,:),:,:),dens(id_s(i_s,:),:,:));
bar(b,h,'FaceAlpha',0.6, 'DisplayName',['M=' num2str(Ma,3) ' ; S=' num2str(Sa,3)]);
legend; xlabel([id_s_name{i_s} '. Total nb: ' num2str(nansum(nansum(nansum(dens(id_s(i_s,:),:,:))))/1000000,3) ' M']);xlim([0 30]);
figure('position',[0 0 400 1000]);
[b,h,Mg,Sg] = histvdens(v_g(id_sp,:,:),dens(id_sp,:,:),0:1:30);
bar(b,h,1,'FaceAlpha',0.6, 'DisplayName',['Spring. M=' num2str(Mg,3) ' ; S=' num2str(Sg,3)],'FaceColor',[.1 .1 .1],'EdgeAlpha',0);
[b,h,Mg,Sg] = histvdens(v_g(id_au,:,:),dens(id_au,:,:),0:1:30);
bar(b,h,1,'FaceAlpha',0.6, 'DisplayName',['Autumn. M=' num2str(Mg,3) ' ; S=' num2str(Sg,3)],'FaceColor',[.5 .5 .5],'EdgeAlpha',0);
legend; xlabel(['Groundspeed [m/s]. Total nb: ' num2str(nansum(nansum(nansum(dens(id_sp,:,:))))/1000000,3) ' M']);
ylabel('Sum of bird density [bird/km^3]'); xlim([0 30]);
[b,h,Ma,Sa] = histvdens(sd_vvp(id_sp,:,:),dens(id_sp,:,:),1.5:0.2:6);
bar(b,h,1,'FaceAlpha',0.6, 'DisplayName',['Spring. M=' num2str(Ma,3) ' ; S=' num2str(Sa,3)],'FaceColor',[.1 .1 .1],'EdgeAlpha',0);
[b,h,Ma,Sa] = histvdens(sd_vvp(id_au,:,:),dens(id_au,:,:),1.5:0.2:6);
bar(b,h,1,'FaceAlpha',0.6, 'DisplayName',['Autumn. M=' num2str(Ma,3) ' ; S=' num2str(Sa,3)],'FaceColor',[.5 .5 .5],'EdgeAlpha',0);
legend; xlabel(['Spectrum width [m/s]. Total nb: ' num2str(nansum(nansum(nansum(dens(id_au,:,:))))/1000000,3) ' M']);
ylabel('Sum of bird density [bird/km^3]'); xlim([1.5 6]);
[b,h,Ma,Sa] = histvdens(v_a(id_sp,:,:),dens(id_sp,:,:),5:1:20);
bar(b,h,1,'FaceAlpha',0.6, 'DisplayName',['Spring. M=' num2str(Ma,3) ' ; S=' num2str(Sa,3)],'FaceColor',[.1 .1 .1],'EdgeAlpha',0);
[b,h,Ma,Sa] = histvdens(v_a(id_au,:,:),dens(id_au,:,:),5:1:20);
bar(b,h,1,'FaceAlpha',0.6, 'DisplayName',['Autumn. M=' num2str(Ma,3) ' ; S=' num2str(Sa,3)],'FaceColor',[.5 .5 .5],'EdgeAlpha',0);
legend; xlabel(['Airspeed [m/s]. Total nb: ' num2str(nansum(nansum(nansum(dens(id_au,:,:))))/1000000,3) ' M']);
ylabel('Sum of bird density [bird/km^3]'); xlim([5 20]);
figure('position',[0 0 1200 600]);
[b,h,Mg,Sg] = histvdens(v_g(id_sp,:,id_so),dens(id_sp,:,id_so));
bar(b,h,'FaceAlpha',0.6, 'DisplayName', ['Groundspeed. M=' num2str(Mg,3) ' ; S=' num2str(Sg,3)]);
[b,h,Ma,Sa] = histvdens(v_a(id_sp,:,id_so),dens(id_sp,:,id_so));
bar(b,h,'FaceAlpha',0.6, 'DisplayName', ['Airspeed. M=' num2str(Ma,3) ' ; S=' num2str(Sa,3)]);
[b,h,Mw,Sw] = histvdens(v_w(id_sp,:,id_so),dens(id_sp,:,id_so));
yyaxis right; bar(b,h,'FaceAlpha',0.6, 'DisplayName', ['Windspeed. M=' num2str(Mw,3) ' ; S=' num2str(Sw,3)]);
[b,h,Ms,Ss] = histvdens(sd_vvp(id_sp,:,id_so),dens(id_sp,:,id_so));
yyaxis left; bar(b,h,'FaceAlpha',0.6, 'DisplayName', ['sd_{vvp}. M=' num2str(Ms,3) ' ; S=' num2str(Ss,3)]);
ax = gca; ax.YAxis(1).Color = 'k'; ax.YAxis(2).Color = 'k';
legend; xlabel(['South-Spring. Total nb: ' num2str(nansum(nansum(nansum(dens(id_sp,:,id_so),3)))/1000000,3) ' M']);
ylabel('Sum of bird density [bird/km^3]'); xlim([0 30]);
[b,h,Mg,Sg] = histvdens(v_g(id_au,:,id_so),dens(id_au,:,id_so));
bar(b,h,'FaceAlpha',0.6, 'DisplayName', ['Groundspeed. M=' num2str(Mg,3) ' ; S=' num2str(Sg,3)]);
[b,h,Ma,Sa] = histvdens(v_a(id_au,:,id_so),dens(id_au,:,id_so));
bar(b,h,'FaceAlpha',0.6, 'DisplayName', ['Airspeed. M=' num2str(Ma,3) ' ; S=' num2str(Sa,3)]);
[b,h,Mw,Sw] = histvdens(v_w(id_au,:,id_so),dens(id_au,:,id_so));
yyaxis right; bar(b,h,'FaceAlpha',0.6, 'DisplayName', ['Windspeed. M=' num2str(Mw,3) ' ; S=' num2str(Sw,3)]);
[b,h,Ms,Ss] = histvdens(sd_vvp(id_au,:,id_so),dens(id_au,:,id_so));
yyaxis left; bar(b,h,'FaceAlpha',0.6, 'DisplayName', ['sd_{vvp}. M=' num2str(Ms,3) ' ; S=' num2str(Ss,3)]);
ax = gca; ax.YAxis(1).Color = 'k'; ax.YAxis(2).Color = 'k';xlim([0 30])
legend; xlabel(['South-Autumn. Total nb: ' num2str(nansum(nansum(nansum(dens(id_au,:,id_so),3)))/1000000,3) ' M']);
[b,h,Mg,Sg] = histvdens(v_g(id_sp,:,id_no),dens(id_sp,:,id_no));
bar(b,h,'FaceAlpha',0.6, 'DisplayName', ['Groundspeed. M=' num2str(Mg,3) ' ; S=' num2str(Sg,3)]);
[b,h,Ma,Sa] = histvdens(v_a(id_sp,:,id_no),dens(id_sp,:,id_no));
bar(b,h,'FaceAlpha',0.6, 'DisplayName', ['Airspeed. M=' num2str(Ma,3) ' ; S=' num2str(Sa,3)]);
[b,h,Mw,Sw] = histvdens(v_w(id_sp,:,id_no),dens(id_sp,:,id_no));
yyaxis right; bar(b,h,'FaceAlpha',0.6, 'DisplayName', ['Windspeed. M=' num2str(Mw,3) ' ; S=' num2str(Sw,3)]);
[b,h,Ms,Ss] = histvdens(sd_vvp(id_sp,:,id_no),dens(id_sp,:,id_no));
yyaxis left; bar(b,h,'FaceAlpha',0.6, 'DisplayName', ['sd_{vvp}. M=' num2str(Ms,3) ' ; S=' num2str(Ss,3)]);
ax = gca; ax.YAxis(1).Color = 'k'; ax.YAxis(2).Color = 'k';xlim([0 30])
legend; xlabel(['North-Spring. Total nb: ' num2str(nansum(nansum(nansum(dens(id_sp,:,id_no),3)))/1000000,3) ' M']);
[b,h,Mg,Sg] = histvdens(v_g(id_au,:,id_no),dens(id_au,:,id_no));
bar(b,h,'FaceAlpha',0.6, 'DisplayName', ['Groundspeed. M=' num2str(Mg,3) ' ; S=' num2str(Sg,3)]);
[b,h,Ma,Sa] = histvdens(v_a(id_au,:,id_no),dens(id_au,:,id_no));
bar(b,h,'FaceAlpha',0.6, 'DisplayName', ['Airspeed. M=' num2str(Ma,3) ' ; S=' num2str(Sa,3)]);
[b,h,Mw,Sw] = histvdens(v_w(id_au,:,id_no),dens(id_au,:,id_no));
yyaxis right; bar(b,h,'FaceAlpha',0.6, 'DisplayName', ['Windspeed. M=' num2str(Mw,3) ' ; S=' num2str(Sw,3)]);
[b,h,Ms,Ss] = histvdens(sd_vvp(id_au,:,id_no),dens(id_au,:,id_no));
yyaxis left; bar(b,h,'FaceAlpha',0.6, 'DisplayName', ['sd_{vvp}. M=' num2str(Ms,3) ' ; S=' num2str(Ss,3)]);
ax = gca; ax.YAxis(1).Color = 'k'; ax.YAxis(2).Color = 'k'; xlim([0 30])
legend; xlabel(['North-Autumn. Total nb: ' num2str(nansum(nansum(nansum(dens(id_au,:,id_no),3)))/1000000,3) ' M']);
figure('position',[0 0 1600 600]);
subplot(2,4,i_s); hold on;
[b,h,Mg,Sg] = histvdens(v_g(id_s(i_s,:),:,id_fr),dens(id_s(i_s,:),:,id_fr));
bar(b,h,'FaceAlpha',0.6);
[b,h,Ma,Sa] = histvdens(v_a(id_s(i_s,:),:,id_fr),dens(id_s(i_s,:),:,id_fr));
bar(b,h,'FaceAlpha',0.6);
if i_s==1, ylabel('France'), end
[b,h,Mw,Sw] = histvdens(v_w(id_s(i_s,:),:,id_fr),dens(id_s(i_s,:),:,id_fr));
ax = gca; set(ax,'yticklabel',[]); yyaxis right; bar(b,h,'FaceAlpha',0.6);
ax = gca; set(ax,'yticklabel',[]); ax.YAxis(1).Color = 'k'; ax.YAxis(2).Color = 'k';
legend(['M=' num2str(Mg,3) ' ; S=' num2str(Sg,3)], ['M=' num2str(Ma,3) ' ; S=' num2str(Sa,3)], ['M=' num2str(Mw,3) ' ; S=' num2str(Sw,3)])
xlabel(['Total nb: ' num2str(nansum(nansum(nansum(dens(id_s(i_s,:),:,id_fr),3)))/1000000,3) ' M']);xlim([0 30]); box on;
subplot(2,4,i_s+4); hold on;
[b,h,Mg,Sg] = histvdens(v_g(id_s(i_s,:),:,id_de),dens(id_s(i_s,:),:,id_de));
bar(b,h,'FaceAlpha',0.6);
[b,h,Ma,Sa] = histvdens(v_a(id_s(i_s,:),:,id_de),dens(id_s(i_s,:),:,id_de));
bar(b,h,'FaceAlpha',0.6);
[b,h,Mw,Sw] = histvdens(v_w(id_s(i_s,:),:,id_de),dens(id_s(i_s,:),:,id_de));
if i_s==1, ylabel('Germany'), end
ax = gca; set(ax,'yticklabel',[]); yyaxis right; bar(b,h,'FaceAlpha',0.6);
ax = gca; set(ax,'yticklabel',[]); ax.YAxis(1).Color = 'k'; ax.YAxis(2).Color = 'k';
legend(['M=' num2str(Mg,3) ' ; S=' num2str(Sg,3)], ['M=' num2str(Ma,3) ' ; S=' num2str(Sa,3)], ['M=' num2str(Mw,3) ' ; S=' num2str(Sw,3)])
xlabel([id_s_name{i_s} '. Total nb: ' num2str(nansum(nansum(nansum(dens(id_s(i_s,:),:,id_de),3)))/1000000,3) ' M']);xlim([0 30]); box on;
Mg=nan(numel(dc),size(id_s,1)); Mgw=Mg; gu=Mg; gv=Mg; gwu=Mg; gwv=Mg; awu=Mg; awv=Mg; wwu=Mg; wwv=Mg; Maw=Mg; Mww=Mg;
Md(i_d,i_s) = nanmean(reshape(dens(id_s(i_s,:),:,i_d),1,[]));
Mg(i_d,i_s) = nanmean(reshape( v_g(id_s(i_s,:),:,i_d) ,1,[]));
gu(i_d,i_s) = nanmean(reshape( real( gs(id_s(i_s,:),:,i_d)),1,[]));
gv(i_d,i_s) = nanmean(reshape( imag( gs(id_s(i_s,:),:,i_d)),1,[]));
Mw(i_d,i_s) = nanmean(reshape( v_w(id_s(i_s,:),:,i_d) ,1,[]));
wu(i_d,i_s) = nanmean(reshape( real( ws(id_s(i_s,:),:,i_d)),1,[]));
wv(i_d,i_s) = nanmean(reshape( imag( ws(id_s(i_s,:),:,i_d)),1,[]));
[~,~,Mgw(i_d,i_s)] = histvdens( v_g(id_s(i_s,:),:,i_d) ,dens(id_s(i_s,:),:,i_d));
[~,~,gwu(i_d,i_s)] = histvdens( real( gs(id_s(i_s,:),:,i_d)),dens(id_s(i_s,:),:,i_d));
[~,~,gwv(i_d,i_s)] = histvdens( imag( gs(id_s(i_s,:),:,i_d)),dens(id_s(i_s,:),:,i_d));
[~,~,Maw(i_d,i_s)] = histvdens( v_a(id_s(i_s,:),:,i_d) ,dens(id_s(i_s,:),:,i_d));
[~,~,awu(i_d,i_s)] = histvdens( real( as(id_s(i_s,:),:,i_d)),dens(id_s(i_s,:),:,i_d));
[~,~,awv(i_d,i_s)] = histvdens( imag( as(id_s(i_s,:),:,i_d)),dens(id_s(i_s,:),:,i_d));
[~,~,Mww(i_d,i_s)] = histvdens( v_w(id_s(i_s,:),:,i_d) ,dens(id_s(i_s,:),:,i_d));
[~,~,wwu(i_d,i_s)] = histvdens( real( ws(id_s(i_s,:),:,i_d)),dens(id_s(i_s,:),:,i_d));
[~,~,wwv(i_d,i_s)] = histvdens( imag( ws(id_s(i_s,:),:,i_d)),dens(id_s(i_s,:),:,i_d));
figure('position',[0 0 1600 2400]);
subplot(6,4,i_s); hold on;
h=worldmap([floor(min([dc.lat])) ceil(max([dc.lat]))], [floor(min([dc.lon])) ceil(max([dc.lon]))]);
setm(h,'frame','on','grid','off'); set(findall(h,'Tag','MLabel'),'visible','off'); set(findall(h,'Tag','PLabel'),'visible','off')
plotm(coastlat, coastlon,'k'); % geoshow('worldrivers.shp','Color', 'blue');
scatterm([dc(id_de).lat]',[dc(id_de).lon]',100,Md(id_de,i_s),'filled');
scatterm([dc(~id_de&~id_fr).lat]',[dc(~id_de&~id_fr).lon]',100,Md(~id_de&~id_fr,i_s),'s','filled');
scatterm([dc(id_fr).lat]',[dc(id_fr).lon]',100,Md(id_fr,i_s),'d','filled');
if i_s==1, title(['Density - ' id_s_name{i_s}]); end
subplot(6,4,i_s+4); hold on;
h=worldmap([floor(min([dc.lat])) ceil(max([dc.lat]))], [floor(min([dc.lon])) ceil(max([dc.lon]))]);
setm(h,'frame','on','grid','off'); set(findall(h,'Tag','MLabel'),'visible','off'); set(findall(h,'Tag','PLabel'),'visible','off')
plotm(coastlat, coastlon,'k'); % geoshow('worldrivers.shp','Color', 'blue');
scatterm([dc(id_de).lat]',[dc(id_de).lon]',100,Mw(id_de,i_s),'filled');
scatterm([dc(~id_de&~id_fr).lat]',[dc(~id_de&~id_fr).lon]',100,Mw(~id_de&~id_fr,i_s),'s','filled');
scatterm([dc(id_fr).lat]',[dc(id_fr).lon]',100,Mw(id_fr,i_s),'d','filled');
quiverm([dc.lat]',[dc.lon]',wu(:,i_s),wv(:,i_s),'k')
if i_s==1, title('All windspeed'); end
subplot(6,4,i_s+8); hold on;
h=worldmap([floor(min([dc.lat])) ceil(max([dc.lat]))], [floor(min([dc.lon])) ceil(max([dc.lon]))]);
setm(h,'frame','on','grid','off'); set(findall(h,'Tag','MLabel'),'visible','off'); set(findall(h,'Tag','PLabel'),'visible','off')
plotm(coastlat, coastlon,'k'); % geoshow('worldrivers.shp','Color', 'blue');
scatterm([dc(id_de).lat]',[dc(id_de).lon]',100,Mg(id_de,i_s),'filled');
scatterm([dc(~id_de&~id_fr).lat]',[dc(~id_de&~id_fr).lon]',100,Mg(~id_de&~id_fr,i_s),'s','filled');
scatterm([dc(id_fr).lat]',[dc(id_fr).lon]',100,Mg(id_fr,i_s),'d','filled');
quiverm([dc.lat]',[dc.lon]',gu(:,i_s),gv(:,i_s),'k')
if i_s==1, title('All groundspeed'); end
subplot(6,4,i_s+12); hold on;
h=worldmap([floor(min([dc.lat])) ceil(max([dc.lat]))], [floor(min([dc.lon])) ceil(max([dc.lon]))]);
setm(h,'frame','on','grid','off'); set(findall(h,'Tag','MLabel'),'visible','off'); set(findall(h,'Tag','PLabel'),'visible','off')
plotm(coastlat, coastlon,'k'); % geoshow('worldrivers.shp','Color', 'blue');
scatterm([dc(id_de).lat]',[dc(id_de).lon]',100,Mww(id_de,i_s),'filled');
scatterm([dc(~id_de&~id_fr).lat]',[dc(~id_de&~id_fr).lon]',100,Mww(~id_de&~id_fr,i_s),'s','filled');
scatterm([dc(id_fr).lat]',[dc(id_fr).lon]',100,Mww(id_fr,i_s),'d','filled');
quiverm([dc.lat]',[dc.lon]',wwu(:,i_s),wwv(:,i_s),'k')
if i_s==1, title('Weighted windspeed'); end
subplot(6,4,i_s+16); hold on;
h=worldmap([floor(min([dc.lat])) ceil(max([dc.lat]))], [floor(min([dc.lon])) ceil(max([dc.lon]))]);
setm(h,'frame','on','grid','off'); set(findall(h,'Tag','MLabel'),'visible','off'); set(findall(h,'Tag','PLabel'),'visible','off')
plotm(coastlat, coastlon,'k'); % geoshow('worldrivers.shp','Color', 'blue');
scatterm([dc(id_de).lat]',[dc(id_de).lon]',100,Mgw(id_de,i_s),'filled');
scatterm([dc(~id_de&~id_fr).lat]',[dc(~id_de&~id_fr).lon]',100,Mgw(~id_de&~id_fr,i_s),'s','filled');
scatterm([dc(id_fr).lat]',[dc(id_fr).lon]',100,Mgw(id_fr,i_s),'d','filled');
quiverm([dc.lat]',[dc.lon]',gwu(:,i_s),gwv(:,i_s),'k')
if i_s==1, title('Weighted groundspeed'); end
subplot(6,4,i_s+20); hold on;
h=worldmap([floor(min([dc.lat])) ceil(max([dc.lat]))], [floor(min([dc.lon])) ceil(max([dc.lon]))]);
setm(h,'frame','on','grid','off'); set(findall(h,'Tag','MLabel'),'visible','off'); set(findall(h,'Tag','PLabel'),'visible','off')
plotm(coastlat, coastlon,'k'); % geoshow('worldrivers.shp','Color', 'blue');
scatterm([dc(id_de).lat]',[dc(id_de).lon]',100,Maw(id_de,i_s),'filled');
scatterm([dc(~id_de&~id_fr).lat]',[dc(~id_de&~id_fr).lon]',100,Maw(~id_de&~id_fr,i_s),'s','filled');
scatterm([dc(id_fr).lat]',[dc(id_fr).lon]',100,Maw(id_fr,i_s),'d','filled');
quiverm([dc.lat]',[dc.lon]',awu(:,i_s),awv(:,i_s),'k')
if i_s==1, title('Weighted airspeed'); end
figure('position',[0 0 1000 1000]);
eb(1) = errorbar(Mg,Ma,Sg/10, 'horizontal', 'LineStyle', 'none','Color','k');
eb(2) = errorbar(Mg,Ma,Sa/10, 'vertical', 'LineStyle', 'none','Color','k');
scatter(Mg(id_de),Ma(id_de),100,[dc(id_de).lat],'filled');
scatter(Mg(~id_de&~id_fr),Ma(~id_de&~id_fr),100,[dc(~id_de&~id_fr).lat],'s','filled');
scatter(Mg(id_fr),Ma(id_fr),100,[dc(id_fr).lat],'d','filled');
plot([7 10.5],[7 10.5],'-r')
text(Mg,Ma-0.1,{dc.name},'HorizontalAlignment','center')
axis equal; xlabel('Ground Speed');ylabel('Air Speed');
c=colorbar(); c.Label.String = 'Lattitude';
title('Used=average weighted by density')
Mg2=nanmean(reshape(v_g,[],numel(dc)));
Ma2=nanmean(reshape(v_a,[],numel(dc)));
Sg2=nanstd(reshape(v_g,[],numel(dc)));
Sa2=nanstd(reshape(v_a,[],numel(dc)));
eb(1) = errorbar(Mg2,Ma2,Sg2/10, 'horizontal', 'LineStyle', 'none','Color','k');
eb(2) = errorbar(Mg2,Ma2,Sa2/10, 'vertical', 'LineStyle', 'none','Color','k');
scatter(Mg2(id_de),Ma2(id_de),100,[dc(id_de).lat],'filled');
scatter(Mg2(~id_de&~id_fr),Ma2(~id_de&~id_fr),100,[dc(~id_de&~id_fr).lat],'s','filled');
scatter(Mg2(id_fr),Ma2(id_fr),100,[dc(id_fr).lat],'d','filled');
plot([7 10.5],[7 10.5],'-r')
text(Mg2,Ma2-0.1,{dc.name},'HorizontalAlignment','center')
axis equal; xlabel('Ground Speed');ylabel('Air Speed');
c=colorbar(); c.Label.String = 'Lattitude'; title('Available')
a = [0 600 1200 2000 5000];
figure('position',[0 0 900 600]);
id_a = dc(1).alt>a(i_a) & dc(1).alt<=a(i_a+1);
subplot(numel(a)-1,1,numel(a)-i_a); hold on
[b,h,Mg,Sg] = histvdens(v_g(id_sp,id_a,:),dens(id_sp,id_a,:));
bar(b,h,'FaceAlpha',0.6);
[b,h,Ma,Sa] = histvdens(v_a(id_sp,id_a,:),dens(id_sp,id_a,:));
bar(b,h,'FaceAlpha',0.6);
[b,h,Mw,Sw] = histvdens(v_w(id_sp,id_a,:),dens(id_sp,id_a,:));
yyaxis right; bar(b,h,'FaceAlpha',0.6);
[b,h, Ms, Ss ] = histvdens(sd_vvp(id_sp,id_a,:),dens(id_sp,id_a,:));
yyaxis left; bar(b,h,'FaceAlpha',0.6);
ax = gca; ax.YAxis(1).Color = 'k'; ax.YAxis(2).Color = 'k';
legend(['Groundspeed. M=' num2str(Mg,3) ' ; S=' num2str(Sg,3)], ['Airspeed. M=' num2str(Ma,3) ' ; S=' num2str(Sa,3)], ...
['SD VVP. M=' num2str(Ms,3) ' ; S=' num2str(Ss,3)], ['Windspeed. M=' num2str(Mw,3) ' ; S=' num2str(Sw,3)])
xlabel(['Spring. Alt:' num2str(min(dc(1).alt(id_a))) 'm - ' num2str(max(dc(1).alt(id_a))) 'm. Total: ' num2str(nansum(nansum(nansum(dens(id_sp,id_a,:),3)))/1000000,2) ' M']); xlim([0 30]);
a = [0 600 1200 2000 5000];
figure('position',[0 0 900 600]);
id_a = dc(1).alt>a(i_a) & dc(1).alt<=a(i_a+1);
subplot(numel(a)-1,1,numel(a)-i_a); hold on
[b,h,Mg,Sg] = histvdens(v_g(id_au,id_a,:),dens(id_au,id_a,:));
bar(b,h,'FaceAlpha',0.6);
[b,h,Ma,Sa] = histvdens(v_a(id_au,id_a,:),dens(id_au,id_a,:));
bar(b,h,'FaceAlpha',0.6);
[b,h,Mw,Sw] = histvdens(v_w(id_au,id_a,:),dens(id_au,id_a,:));
yyaxis right; bar(b,h,'FaceAlpha',0.6);
[b,h, Ms, Ss] = histvdens(sd_vvp(id_au,id_a,:),dens(id_au,id_a,:));
yyaxis left; bar(b,h,'FaceAlpha',0.6);
ax = gca; ax.YAxis(1).Color = 'k'; ax.YAxis(2).Color = 'k';
legend(['Groundspeed. M=' num2str(Mg,3) ' ; S=' num2str(Sg,3)], ['Airspeed. M=' num2str(Ma,3) ' ; S=' num2str(Sa,3)], ...
['SD VVP. M=' num2str(Ms,3) ' ; S=' num2str(Ss,3)], ['Windspeed. M=' num2str(Mw,3) ' ; S=' num2str(Sw,3)])
xlabel(['Autumn. Alt:' num2str(min(dc(1).alt(id_a))) 'm - ' num2str(max(dc(1).alt(id_a))) 'm. Total: ' num2str(nansum(nansum(nansum(dens(id_au,id_a,:),3)))/1000000,2) ' M']); xlim([0 30]);
Mg = nan(numel(dc(1).alt),2);Sg=Mg; Ma=Mg; Sa=Mg; Mw=Mg; Sw=Mg; Ms=Mg; Ss=Mg; Md = Mg; Sd=Mg;
for i_a= 1:numel(dc(1).alt)
[~,~,Mg(i_a,1), Sg(i_a,1)] = histvdens(v_g(id_sp,i_a,:),dens(id_sp,i_a,:));
[~,~,Ma(i_a,1), Sa(i_a,1)] = histvdens(v_a(id_sp,i_a,:),dens(id_sp,i_a,:));
[~,~,Mw(i_a,1), Sw(i_a,1)] = histvdens(v_w(id_sp,i_a,:),dens(id_sp,i_a,:));
[~,~,Ms(i_a,1), Ss(i_a,1)] = histvdens(sd_vvp(id_sp,i_a,:),dens(id_sp,i_a,:));
[~,~,Md(i_a,1), Sd(i_a,1)] = histvdens(dens(id_sp,i_a,:),dens(id_sp,i_a,:));
[~,~,Mg(i_a,2), Sg(i_a,2)] = histvdens(v_g(id_au,i_a,:),dens(id_au,i_a,:));
[~,~,Ma(i_a,2), Sa(i_a,2)] = histvdens(v_a(id_au,i_a,:),dens(id_au,i_a,:));
[~,~,Mw(i_a,2), Sw(i_a,2)] = histvdens(v_w(id_au,i_a,:),dens(id_au,i_a,:));
[~,~,Ms(i_a,2), Ss(i_a,2)] = histvdens(sd_vvp(id_au,i_a,:),dens(id_au,i_a,:));
[~,~,Md(i_a,2), Sd(i_a,2)] = histvdens(dens(id_au,i_a,:),dens(id_au,i_a,:));
figure('position',[0 0 900 600]);
subplot(1,3,[1 2]); hold on;
tmp = errorbar(Mg(:,1),dc(1).alt,Sg(:,1)/10, 'horizontal', 'LineStyle', 'none'); l(1)=plot(Mg(:,1),dc(1).alt,'-','Color',tmp.Color,'LineWidth',2);
tmp = errorbar(Ma(:,1),dc(1).alt,Sa(:,1)/10, 'horizontal', 'LineStyle', 'none'); l(2)=plot(Ma(:,1),dc(1).alt,'-','Color',tmp.Color,'LineWidth',2);
tmp = errorbar(Mw(:,1),dc(1).alt,Sw(:,1)/10, 'horizontal', 'LineStyle', 'none'); l(3)=plot(Mw(:,1),dc(1).alt,'-','Color',tmp.Color,'LineWidth',2);
tmp = errorbar(Ms(:,1),dc(1).alt,Ss(:,1)/10, 'horizontal', 'LineStyle', 'none'); l(4)=plot(Ms(:,1),dc(1).alt,'-','Color',tmp.Color,'LineWidth',2);
xlim([0 20]);ylabel('Altitude asl [m]'); xlabel('Spring. Speed [m/s]'); ylim([0 4000]);
ax1=gca; ax2 = axes('Position',ax1.Position,'XAxisLocation','top','Color','none','YTickLabel',[]); hold on; c = get(gca,'colororder');
tmp = errorbar(Md(:,1),dc(1).alt,Sd(:,1)/10, 'horizontal', 'LineStyle', 'none','Parent', ax2,'Color',c(6,:)); l(5)=plot(Md(:,1),dc(1).alt,'-','Color',tmp.Color,'LineWidth',2,'Parent',ax2,'Color',c(6,:));
xlabel('Density [bird/km^3]'); ylim([0 4000]); xlim([0 120]); legend(l,{'Groundspeed','Airspeed','Windspeed','SD VVP','Density'});
tmp = errorbar(Mg(:,2),dc(1).alt,Sg(:,2)/10, 'horizontal', 'LineStyle', 'none'); l(1)=plot(Mg(:,2),dc(1).alt,'-','Color',tmp.Color,'LineWidth',2);
tmp = errorbar(Ma(:,2),dc(1).alt,Sa(:,2)/10, 'horizontal', 'LineStyle', 'none'); l(2)=plot(Ma(:,2),dc(1).alt,'-','Color',tmp.Color,'LineWidth',2);
tmp = errorbar(Mw(:,2),dc(1).alt,Sw(:,2)/10, 'horizontal', 'LineStyle', 'none'); l(3)=plot(Mw(:,2),dc(1).alt,'-','Color',tmp.Color,'LineWidth',2);
tmp = errorbar(Ms(:,2),dc(1).alt,Ss(:,2)/10, 'horizontal', 'LineStyle', 'none'); l(4)=plot(Ms(:,2),dc(1).alt,'-','Color',tmp.Color,'LineWidth',2);
xlim([3 12]); xlabel('Autumn. Speed [m/s]'); ylim([0 4000]); ax1=gca; ax2 = axes('Position',ax1.Position,'XAxisLocation','top','Color','none','YTickLabel',[]); hold on; c = get(gca,'colororder');
tmp = errorbar(Md(:,2),dc(1).alt,Sd(:,2)/10, 'horizontal', 'LineStyle', 'none','Parent', ax2,'Color',c(6,:)); l(5)=plot(Md(:,2),dc(1).alt,'-','Color',tmp.Color,'LineWidth',2,'Parent',ax2,'Color',c(6,:));
xlabel('Density [bird/km^3]'); ylim([0 4000]); xlim([0 60])
Mg = nan(numel(dc(1).alt),4);Sg=Mg; Ma=Mg; Sa=Mg; Mw=Mg; Sw=Mg; Ms=Mg; Ss=Mg; Md = Mg; Sd=Mg;
for i_a= 1:numel(dc(1).alt)
[~,~,Mg(i_a,i_s), Sg(i_a,i_s)] = histvdens(v_g(id_s(i_s,:),i_a,:),dens(id_s(i_s,:),i_a,:));
[~,~,Ma(i_a,i_s), Sa(i_a,i_s)] = histvdens(v_a(id_s(i_s,:),i_a,:),dens(id_s(i_s,:),i_a,:));
[~,~,Mw(i_a,i_s), Sw(i_a,i_s)] = histvdens(v_w(id_s(i_s,:),i_a,:),dens(id_s(i_s,:),i_a,:));
[~,~,Ms(i_a,i_s), Ss(i_a,i_s)] = histvdens(sd_vvp(id_s(i_s,:),i_a,:),dens(id_s(i_s,:),i_a,:));
[~,~,Md(i_a,i_s), Sd(i_a,i_s)] = histvdens(dens(id_s(i_s,:),i_a,:),dens(id_s(i_s,:),i_a,:));
figure('position',[0 0 1600 600]);
subplot(1,4,i_s); hold on;
tmp = errorbar(Mg(:,i_s),dc(i_s).alt,Sg(:,i_s)/10, 'horizontal', 'LineStyle', 'none'); l(1)=plot(Mg(:,i_s),dc(i_s).alt,'-','Color',tmp.Color,'LineWidth',2);
tmp = errorbar(Ma(:,i_s),dc(i_s).alt,Sa(:,i_s)/10, 'horizontal', 'LineStyle', 'none'); l(2)=plot(Ma(:,i_s),dc(i_s).alt,'-','Color',tmp.Color,'LineWidth',2);
tmp = errorbar(Mw(:,i_s),dc(i_s).alt,Sw(:,i_s)/10, 'horizontal', 'LineStyle', 'none'); l(3)=plot(Mw(:,i_s),dc(i_s).alt,'-','Color',tmp.Color,'LineWidth',2);
tmp = errorbar(Ms(:,i_s),dc(i_s).alt,Ss(:,i_s)/10, 'horizontal', 'LineStyle', 'none'); l(4)=plot(Ms(:,i_s),dc(i_s).alt,'-','Color',tmp.Color,'LineWidth',2);
xlim([0 20]);ylabel('Altitude asl [m]'); xlabel([id_s_name{i_s} '. Speed [m/s]']); ylim([0 4000]);
ax1=gca; ax2 = axes('Position',ax1.Position,'XAxisLocation','top','Color','none','YTickLabel',[]); hold on; c = get(gca,'colororder');
tmp = errorbar(Md(:,i_s),dc(1).alt,Sd(:,i_s)/10, 'horizontal', 'LineStyle', 'none','Parent', ax2,'Color',c(6,:)); l(5)=plot(Md(:,i_s),dc(1).alt,'-','Color',tmp.Color,'LineWidth',2,'Parent',ax2,'Color',c(6,:));
xlabel('Density [bird/km^3]'); ylim([0 4000]); xlim([0 70]); legend(l,{'Groundspeed','Airspeed','Windspeed','SD VVP','Density'});
mida=(a(2:end)+a(1:end-1))/2;
Mg = nan(numel(a)-1,2);Sg=Mg; Ma=Mg; Sa=Mg; Mw=Mg; Sw=Mg; Ms=Mg; Ss=Mg; Md = Mg; Sd=Mg;
id_a = repmat(permute(NNT>a(i_a) & NNT<=a(i_a+1) & repmat(id_sp',1,numel(dc)),[1 3 2]),1,25,1);
[~,~,Mg(i_a,1), Sg(i_a,1)] = histvdens(v_g(id_a),dens(id_a));
[~,~,Ma(i_a,1), Sa(i_a,1)] = histvdens(v_a(id_a),dens(id_a));
[~,~,Mw(i_a,1), Sw(i_a,1)] = histvdens(v_w(id_a),dens(id_a));
[~,~,Ms(i_a,1), Ss(i_a,1)] = histvdens(sd_vvp(id_a),dens(id_a));
[~,~,Md(i_a,1), Sd(i_a,1)] = histvdens(dens(id_a),dens(id_a));
id_a = repmat(permute(NNT>a(i_a) & NNT<=a(i_a+1) & repmat(id_au',1,numel(dc)),[1 3 2]),1,25,1);
[~,~,Mg(i_a,2), Sg(i_a,2)] = histvdens(v_g(id_a),dens(id_a));
[~,~,Ma(i_a,2), Sa(i_a,2)] = histvdens(v_a(id_a),dens(id_a));
[~,~,Mw(i_a,2), Sw(i_a,2)] = histvdens(v_w(id_a),dens(id_a));
[~,~,Ms(i_a,2), Ss(i_a,2)] = histvdens(sd_vvp(id_a),dens(id_a));
[~,~,Md(i_a,2), Sd(i_a,2)] = histvdens(dens(id_a),dens(id_a));
figure('position',[0 0 900 800]);
tmp = errorbar(mida, Mg(:,1),Sg(:,1)/10, 'vertical', 'LineStyle', 'none'); l(1)=plot(mida, Mg(:,1),'-','Color',tmp.Color,'LineWidth',2);
tmp = errorbar(mida, Ma(:,1),Sa(:,1)/10, 'vertical', 'LineStyle', 'none'); l(2)=plot(mida, Ma(:,1),'-','Color',tmp.Color,'LineWidth',2);
tmp = errorbar(mida, Mw(:,1),Sw(:,1)/10, 'vertical', 'LineStyle', 'none'); l(3)=plot(mida, Mw(:,1),'-','Color',tmp.Color,'LineWidth',2);
tmp = errorbar(mida, Ms(:,1),Ss(:,1)/10, 'vertical', 'LineStyle', 'none'); l(4)=plot(mida, Ms(:,1),'-','Color',tmp.Color,'LineWidth',2);
ylim([3 14]);ylabel('Altitude asl [m]'); ylabel('Spring. Speed [m/s]'); xlim([-1 1]);
ax1=gca; ax2 = axes('Position',ax1.Position,'YAxisLocation','right','Color','none','XTickLabel',[]); hold on; c = get(gca,'colororder');
tmp = errorbar(mida, Md(:,1),Sd(:,1)/10, 'vertical', 'LineStyle', 'none','Parent', ax2,'Color',c(6,:)); l(5)=plot(mida, Md(:,1),'-','Color',tmp.Color,'LineWidth',2,'Parent',ax2,'Color',c(6,:));
ylabel('Density [bird/km^3]'); xlim([-1 1]); ylim([0 50]); legend(l,{'Groundspeed','Airspeed','Windspeed','SD VVP','Density'});
tmp = errorbar(mida, Mg(:,2),Sg(:,2)/10, 'vertical', 'LineStyle', 'none'); l(1)=plot(mida, Mg(:,2),'-','Color',tmp.Color,'LineWidth',2);
tmp = errorbar(mida, Ma(:,2),Sa(:,2)/10, 'vertical', 'LineStyle', 'none'); l(2)=plot(mida, Ma(:,2),'-','Color',tmp.Color,'LineWidth',2);
tmp = errorbar(mida, Mw(:,2),Sw(:,2)/10, 'vertical', 'LineStyle', 'none'); l(3)=plot(mida, Mw(:,2),'-','Color',tmp.Color,'LineWidth',2);
tmp = errorbar(mida, Ms(:,2),Ss(:,2)/10, 'vertical', 'LineStyle', 'none'); l(4)=plot(mida, Ms(:,2),'-','Color',tmp.Color,'LineWidth',2);
ylim([3 14]); ylabel('Autumn. Speed [m/s]'); xlim([-1 1]); xlabel('NNT Progress of Night (-1: Sunset | 1: Sunrise)');
ax1=gca; ax2 = axes('Position',ax1.Position,'YAxisLocation','right','Color','none','XTickLabel',[]); hold on; c = get(gca,'colororder');
tmp = errorbar(mida, Md(:,2),Sd(:,2)/10, 'vertical', 'LineStyle', 'none','Parent', ax2,'Color',c(6,:)); l(5)=plot(mida, Md(:,2),'-','Color',tmp.Color,'LineWidth',2,'Parent',ax2,'Color',c(6,:));
ylabel('Density [bird/km^3]'); xlim([-1 1]); ylim([0 50]);
figure('position',[0 0 900 1400]);
subplot(4,1,i_s); hold on;
tmp = errorbar(mida, Mg(:,i_s),Sg(:,i_s)/10, 'vertical', 'LineStyle', 'none'); l(1)=plot(mida, Mg(:,i_s),'-','Color',tmp.Color,'LineWidth',2);
tmp = errorbar(mida, Ma(:,i_s),Sa(:,i_s)/10, 'vertical', 'LineStyle', 'none'); l(2)=plot(mida, Ma(:,i_s),'-','Color',tmp.Color,'LineWidth',2);
tmp = errorbar(mida, Mw(:,i_s),Sw(:,i_s)/10, 'vertical', 'LineStyle', 'none'); l(3)=plot(mida, Mw(:,i_s),'-','Color',tmp.Color,'LineWidth',2);
tmp = errorbar(mida, Ms(:,i_s),Ss(:,i_s)/10, 'vertical', 'LineStyle', 'none'); l(4)=plot(mida, Ms(:,i_s),'-','Color',tmp.Color,'LineWidth',2);
ylim([3 16]);ylabel('Altitude asl [m]'); ylabel([id_s_name{i_s} '. Speed [m/s]']); xlim([-1 1]);
ax1=gca; ax2 = axes('Position',ax1.Position,'YAxisLocation','right','Color','none','XTickLabel',[]); hold on; c = get(gca,'colororder');
tmp = errorbar(mida, Md(:,i_s),Sd(:,i_s)/10, 'vertical', 'LineStyle', 'none','Parent', ax2,'Color',c(6,:)); l(5)=plot(mida, Md(:,i_s),'-','Color',tmp.Color,'LineWidth',2,'Parent',ax2,'Color',c(6,:));
ylabel('Density [bird/km^3]'); xlim([-1 1]); ylim([0 50]); legend(l,{'Groundspeed','Airspeed','Windspeed','SD VVP','Density'});